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Introduction 

Because of the magnitude of Che accomplishments of Che last 6 months, 
it was necessary to deviate from the usual letter-type report. The 
major accomplishments are: 

. Resolution of Che difficulties in solution methodology to be 
used to deal with the potentially highly nonlinear rotor equa- 
tions when dynamic coupling is included. 

. A solution methodology is selected to solve the nonlinear 

differential equations. The selected method was verified to 
give good results even at large nonlinearity levels. 

. Extension of the transfer matrix methodology to the solution 
of nonlinear problems. It is believed that this is a tech- 
nical first. Transfer matrix methods are obstenslbly a 
linear mechanics method. They now have nonlinear capabilities. 
Moreover, these nonlinear analysis capabilities extend into 
the large nonlinear regime. It Is hoped that this work will 
open up new applications of transfer matrices to structural 
mechanics problems. 

These topics will be discussed in the above order In that this 
chronology presents the logical progression to the transfer matrix 
breakthrough. 
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RESOLUTION OF DIFFICULTIES WITH PREVIOUS SOLUTION METHODOLOGY 


The solution technique as reported in the previous report was found 
to be deficient in several areas. Also, it was somewhat erroneously 
termed a harmonic-balance technique, while, in fact, it was more a 
modification of the perturbation technique. The greatest deficiency was 
that it produced accurate results only for cases where the nonlinearity 
was small, like the perturbation technique. Since our motivation for 
dealing with the dynamically coupled rotor problem is that it is sus- 
pected that the nonlinearities are not small. This was a serious short- 
coming. Also, the approximate solution predicted regions of subharmonic 
instability that could not be substantiated by time transient numerical 
solutions. Hence, a superior technique was sought for solving problems 
of this type. 

The technique investigated is indeed the classical harmonic-balance 
technique, but with a slightly different method of actually obtaining 
the solutl 'n. Indeed, the solution of linear ordinary differential 
equations with periodic solutions uses the harmonic-balance technique. 

For a linear ordinary differential equation with a periodic solu- 
tion, one assumes the form of Che solution and substitutes this assumed 
solution into the differential equation. This results in an algebraic 
equation that can be solved for the unknown response amplitude and phase 
relationship. This results in an exact solution for the linear differ- 
ential equation. 
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In dealing with nonlinear ordinary differencial equations, this 
same technique can be used to obtain approximate periodic solutions. 

A form of the solution is assumed and si >stituted into the differ- 
encial equation. Due to the nonlinearicies, this often results in 
an algebraic equation containing trigonometric functions raised Co 
Integer powers. These terms can be reduced to simple harmonic func- 
tions through Che use of trigonometric identities. The resulting 
equation is algebraic in nature, containing trigonometric terms, 
primarily sines and cosines of many frequencies on the right-hand 
side of the equation while containing only terms of a single frequency 
on Che left-hand side. It can be argued that if the resulting equa- 
tion is indeed an Identity, Che coefficients of the sines and cosines 
on both sides of the equation must form and identity. Hence, one equ- 
ates coefficients of like frequencies of sines and cosines, ignoring chose 
Chose harmonic terms not in the assumed solution, and obtains simple 
algebraic expressions which can be solved for the unknown amplitude 
and phase information. Perhaps at this point an example would serve 
to clarify Che above description. 

SOLUTION OF DUFFING *S EQUATION BY THE HARMONIC-BALANCE TECHNIQUE 

Duffing 's equation is representative of a system consisting of a 
mass resting on a spring with a cubic nonlinearity as shown in Fig. 1. 

The equation of motion is 

mx + k(x+€x^) ■ F 8in(wt) . 
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Fig. 1 Harmonic Oscillator with Nonlinear Spring 
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This equation was chosen because it has been thoroughly investigated 
and the results of the hannonic-balance technique could be easily 
verified. The solution methodology is as follows: 

From Newton's Laws, one can argue that one component of the re- 
sponse must be at the excitation frequency, so initially a solu- 
tion of the form 

x(t) ■ x^ sin(ajt) 

is assumed and substituted into the equation of motion. The re- 
sulting equation is 

2 3 3 

-mto sin(ujt) + k(x^ sin(wt) + £x^ sin (cot)) ■ F sin(ojt) . 

The term sin^(cot) is simplified using the trigonometric identity 
as 

ain^(u)t) ■ 1/4(3 sin(cot) - sin(3cot)). 

which yields the following equation of motion 

2 3 

-<i) sin(u)t) + k/m x^^ sin(cot) + ek/4m x^ (3 sin(cot) 

- sin(3ut)) ■ F/m sin(u)t) . 

At this point, one could apply the principle of harmonic balance, 
which would yield 

+ k/m Xj^ + 3k/4m ex^^^ ■ F/m 

This could be solved for the unknown amplitude x^. In doing this, 
one ignores Che sin(3u)C) term. This assumption is the reason the 
solution is called approximate. Since this sln(3u)t) term was Ignored, 
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a better solution can be obtained by, upon the next iteration, let- 
ting the initial assumed solution contain a term of this form. Thus 

one assumes 

x(t) ■ x^ sin(cot) + x^ sin(3ujt) 

and substitute this into the differential equation. Upon substitu- 
tion into the equation of motion thi. yields 
2 

~mu) (Xj^ sin((i)t) + 9x^ sin(3(i)t)) + k(Xj^ sin(oit) + x^ sin(3o)t)) 

3 

+ £k((Xj^ sin(cot) + x^ sin(3oJt)) ) ■ F sin (wt) . 

After simplification and application of the harmonic-balance princi- 
ple, the differential equation is reduced to two simultaneous equa- 
tions in terms of the unknowns x^ and x^. These equations are (letting 

2 2 

y * u) m/k) for coefficients of sin(ut) 

3/4 - 3/4 + 3/2 ~ F/k • 0. 

and for coefficients of sin(3u)t) 

-1/4 + 3/4 Ex^^ + 3/2 £*3 ^x, + (l-gy^lx^ - 0. 

As can be seen, these equations are coupled and nonlinear, hence 
an explicit solution cannot be formed, except in the case of small 
nonlinearities. However, the real roots of these equations repre- 
sent the approximate solution to the differential equation. The 
methodology for finding the roots to these equations will be out- 
lined in a subsequent section. This process could be repeated, 
that is readjusting the assumed form of the solution to account 
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for the harmonic terms that were ignored in the previous step, 
yielding a better answer. However, the algebraic equations be- 
come extremely complex and hence, Che process is usually termin- 
ated after one or two iterations. 

Damping can readily be included in the problem, rendering 
the equation of motion as 

3 

mx + cx + k(x+ex ) * F 3 in(oJt) . 

The only modification Co the solution process is in the 
assumed form of the solution, which in this case is 

x(t) ■ Xj^j^ sin(a>c) + cos(ojt) + x^^^ sin(3(ijC) cosOojt). 

This is preferred over a complex formulation so that all com- 
plex roots of the algebraic equations can be Ignored. As before, 
this assumed form of the solution is substituted into the equation 
of motion, this yielding four algebraic equations In terms of the 
four unknowns *31* *32^' algebra for this 

case Is unwieldy, the computer program FMACUT(*) was used to per- 
form the algebraic manipulations. As before, the final step In 
the solution is to find the appropriate solution to the resulting 
algebraic equations. 

SOLUTION OP THE NONLINEAR ALGEBRAIC EQUATIONS 

In essence, the hamonlc-balaoce technique for obtaining approximate 
solutions to nonlinear ordinary differential equations reduces the solu- 
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cions to finding the real roots of simultaneous, nonlinear algebraic 
equatiois. Once these several roots are found, a physical principle 
must be employed to select the correct one. Since the equations are 
nonlinear, there is no direct solution technique for finding the roots. 
Hence, one must resort to an iterative technique. The procedure used 
is based on an approximate solution and employs a truncated Taylor 
series to converge to a more precise answer. Again, perhaps the pro- 
cedure is best described with the use of an example. 

The algebraic equations resulting from the analysis of Duffing 's 
equation in the previous section were: 

f(Xj^,x^) - 3/4 ex^^ - 3/4 ex^^x^ + 3/2 ex^x^^ + (I-y^)x^ - F/k - 0 

g(x^,X2> - -1/4 ex^^ + 3/4 exj^ + 3/2 + (1-9 y^)x3 - 0. 

Assuming one knows an approximate solution (this can be a very 

crude approxmation) , these equations can be expanded in a Taylor Series 
about this approximate solution 

f(x^,X3> - f(a^,a^) + 3f/3x^ (x^-a^) + af/ax^ (x^-a^) + .... 0 

and 

gCx^.x^) - g(a^,a3) + ag/axj^ (x^-a^) + ag/ax^ (Xj-a^) + ... - 0 

and by ignoring the higher order terns, these reduce to a set of linear 
algebraic equations which can be solved for a more precise value of the 
solution x^,x^. This process is repeated until a desired accuracy is 
reached. 

The approximate solution used to start this procedure can often be 
a very poor approximation. Usually, this simply requires additional 
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Iterations to converge on the precise value. Sometimes, however, the 
process fails to converge. This doesn't present a real problem in that 
one can specify a limiting number of iterations to be performed. If 
convergence is not achieved within this limiting number of iterations, 
a new starting point is chosen and the process repeated. 

Another problem occurs when the equations have mulitple sets of 
real roots. In Duffing 's Equation, these multiple roots represent the 
well known jump phenomena. Since each set of solutions satisfies the 
equations equally as well, a criteria must be established in order 
to choose the appropriate solution. While further investigation is 
required, at this time it appears that the correct solution for our 
physical systems is the solution that minimizes the total system energy. 

COMPARISON OF APPROXIMATE SOLUTION WITH NUMERICAL RESULTS 

In order to check the accuracy of the harmonlc-balance solution to 
Duffing 's equation, numerical simulations were performed using the IBM 
program CSMP. Since the approximate solution doesn't conalder initial 
conditions, the nutaerical simulation was continued for sufficient periods 
of time to ensure that all transient responses had been eliminated. The 
output of the CSMP simulation was used in a Fourier Analysis routine to 
determine the amplitude, phase and frequency characteristics of the steady* 
state response. The results of several analyses, along with the approximate 
solutions are shown as amplitude and phaae plots shown In Fig. 2 through 
Fig. 5. Mote that the approximate solutions are represented by a smooth 
curve, while the numerical solutions are represented by discrete points. 
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Fig. 4 Steady-Stata Solution of Duffing '$ Equation: Amplitude and 
Phaae Rclationahlp of the Fundamental Ilacmonic Component 
for the Case of Large Nonlinearity. 
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This is due to the fact that the numerical solutions where somewhat 
expensive, hence only a few frequencies were chosen at which to per- 
form the numerical simulations. 

The system represented by Fig. 2 and Fig. 3 had the parameters ad- 
justed so that the nonlinearities were small. As can be seen, the approxi- 
mate and numerical solutions are in good agreement for this case. The 
system represented by Fig. 4 and Fig. 5 had the parameters adjusted so 
that the nonlinearities were large. By large, one means that the maximum 
amplitude of the nonlinear terms in the differential equation were of 
the same magnitude as that of the linear terms. As can be seen, the 
harmonic-balance solution produced reasonably good results in this case. 

As a final check of the harmonic-balance technique, the magnitudes 
of the Fourier components obtained from the numerical simulations 
were compared. These Fourier components appear in Table 1. As can be 
seen from Table 1, only the first and third harmonic components are 
significant, the others being much smaller than these two dor lant 
components. Hence, additional harmonic terms were not necessary in 
the initial assumed solution. 

Hence, the results of this portion of the investigation indicate 
that the harmonic-balance technique can be used to produce reasonably 
good solutions to this type of nonlinear differential equation, even 
when the nonlinearities are large. 
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x(t) » T— + I X cos(ntjJt - f ) 
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FOURIER 

COEFFICIENT 

SMALL 

NONLINEARITY 

LARGE 

NONLINEARITY 


^ -6 

-7 

X 

0 

-3.3 X 10 

-1.0 X 10 

^1 

1.162 

7.38 X 10'^ 


^ « -6 

-7 

^2 

2.8 X 10 

1.0 X 10 


-2 

-1 

^3 

3.61 X 10 

2.308 X 10 

^4 

1.67 X 10“^ 

1.0 X 10*^ 

^5 

7.67 X lO"^ 

2.75 X lO'^ 


TABLE 1 Comparison of Harmonic Components Obtained from the 
Numerical Simulations for a 0.5 Hertz Excitation 
Frequency 
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SOLUTION OF DYNAMICALLY -COUPLED NONLINEAR DIFFERENTIAL EQUATIONS 


The motlvaclon for dealing with this problem is chat the rotor equa- 
tions are not only nonlinear in the sense of having terms raised to 
ii’-eger powers, but they also have terms which contain products of 
t ' ate variables. In Che previous section, the ability of the harmonlc- 
bf;lance technique to solve nonlinear equations of a system with a single 
degree of freedom was demonstrated. Since Che rotor equations are for 
a system with six degrees of freedom and contain the aforementioned non- 
li lear terms, one must demonstrate the ability of the solution technique 
to solve such equations adequately. 

Instead of beginning with the rotor equations, a simpler set of equa- 
tions was chosen. This was to allow us to become familiar with solving 
these types of equations before dealing with the more difficult rotor 
equations. To this end, we chose to analyze Che spring-pendulum system 
shown in Fig. 6. This system was chosen because its equations of motion 
contain terms similar to the rotor equations. Using Lagrange's technique, 
the differential equations of motion for this system can be shown to be 

f - g cos(9) + x/m(r-l) ■ 0 
r^6 + 2rr0 + C./m 9 + gr sin(0) ■ T /m sin(u)t) 

w O 

The tenrs of interest are 

• 2 • • 

r6 and 2rr6 . 

The presence of the trigonometric terms sln(6) and cos(6) pose diffi- 
culty to the solution. Since the rotor equations contain no such crlgono- 
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I 

^ metric terms, a truncated Taylor series was used to approximate them. 

The resulting equations are 

r - r9^ C /m r -g(l-hd^) + k/m(r-l) » 0 
R 

! r^9 + 2rr0 + C^/m 0 + gr(0-0^/6) - T /m sin(cot) 

A comparison of the truncated Taylor series along with the appropriate 
trigonometric functions showed that the approximation was very good 
(within 10%) for angles between + 60 degrees. 

As before, since the torque is driving the system at a particular 
frequency, one component of the response will be at this excitation 
frequency. However, the presence of more than single degree of freedom 
can cause problems in the subsequent solution steps. This is because 
the harmonic-balance technique is very sensitive to the assumed form of 
the solution used in the equations. An assumed solution which contains 
erroneous harmonic terms will result in a poor solution. Hence, one must 
carefully consider the assumed form of the solution used. 

As previously mentioned, one component of the response will be at the 
excitation frequency. Hence, one initially assumes 

I 0(t) ■ 0j^j^ sin(o)t) + cos(a)t) 

I Since there is no forcing function in the r(t) equation (the first 

equation), the r(t) response will be produced by gravity and the coup- 

1 *2 

ling term rd . To gain an insight into the r(t) response, the r(t) 
equation can be rearranged as 


r + k/m(r-l) ■ g(l-*j0^) + r0^. 




This looks like the equation of a harmonic oscillator with the 

2 • 2 • 2 . 
driving functions g(1^^0 ) and rS*". The g term and the r9 term will 

produce a constant response and a response at twice the excitation fre- 
quency, Hence, one assumes the r(t) response to be 

r(t) * r + r^. sin(2cjt) + r.« cos(2cjt), 
o 21 22 

As before, these assumed solutions are submitted into the differ- 
ential equations. By ignoring terms not in the original solution and 
applying the principle of harmonic balance, one produces the five non- 
linear polynomials in terms of the five unknowns* The polynomials, 
which were found using FMACUT, are much too long to be presented lere. 
These polynomials were then used in a solution scheme identical to that 
employed in solving Duffing’s equation. 

Again, the results of the analysis are best shown in terms of 
response graphs. In Fig. 7 through Fig. 9 are the amplitude and phase 
angle for the system of Fig. 6, with system parameters as shown on the 
figures. These system parameters were such that the system reponse 
was fairly large. In light of the nonlinearity in the equations, the 
agreement between the harmonic-balance and numerical solutions is 
considered good. 

These results indicate that the harmonic balance method is indeed 
capable of producing good approximate solutions to coupled nonlinear 
differential equations of this form. The major difficulty is choosing 
the appropriate harmonic terms to be used in the solution. This is 
overcome by a careful examination of the equations and verification by 
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numerical analysis. Also, a good solution was produced using only the 
first terms in the assumed solution; i.e., additional iterations on 
the equations to get higher harmonic components were unnecessary. This 
is promising in that perhaps the solution of the rotor equations will 
not greatly expand the already large transfer matrices. 

EXTENSION OF TRANSFER M ATRIX METHODOLOGY TO SOLVE NONLINEAR PROBLEMS 

Traditionally, transfer matrix methodology has been restricted to 
use in linear system analysis. However, when coupled with the harmonic- 
balance technique previously mentioned, transfer matrices can be used 
to solve nonlinear problems in an approximate fashion. The major dif- 
ference is that the solution process becomes iterative for the nonlinear 
problem. To demonstrate the methodology, an example using springs with 
nonlinear stiffnesses and lumped messes is presented. The matrices 
can then be used to model a system represented by Duffing's equation, 
which provides a check on the accuracy of the technique. 

As with conventional transfer matrix methods, one begins by 
developing the matrices for each individual clement. Hence, we begin 
by developing the matrix for a spring having a force-deflection re- 
lationship of the form: 

F • k(x+€x^) 

where x is the net spring deflection. 
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The free-body diagram of the spring is shown in Fig. 14. The inter- 
nal force is assumed to have two frequency components; i.e., 

N(t) • sln(uit) + Nj sin(3ut) . 

This assumption results from the cubic nonlinearity in the spring. 
Using the force-deflection relationship and assuming a response also of 
the form: 

x(t) • sln(a>t) + sln(3uit) 
yields the following relationship: 

3 

N(t) • k(x^sin(ujt) sln(3wt) • e(Xj^ sin(ut) x^ sln(3o)t)) ). 

Simplifying and applying the principle of harmonic balance to this 
equation produces the following simultaneous equations: 

Nj^/k • Xj^+e(3/4 Xj^^«3/4 x^x^2!2 

and 

Nj/k - *3 + e(-l/4 Xj^^+3/2 Xj^^X3+3/4 x^^) , 

To put these equations into transfer-matrix form, let 
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Defining the state vector as x^ | 1 these equations 

can be written in matrix form as 
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where 

■ -e(3/4 Xj^^-3/4 Xj^^x^+3/2 

NLj ■ -e(-l/4 Xj^^3/2 Xj^^X2+3/-+ x^^) 

Note that the nonlinear terms are not separated Into left and right 
variables, but simply carried along in the extension column as a 
correction term. This transfer matrix is best understood as a set 
of nonlinear polynomials to which a solution must be found, in cont ast 
t<> a standard transfer matrix which represents a set of simultaneous 
linear equations which can be explicitly solved. Hence, the solution 
using this matrix in effect seeks values of the system variables 
which satisfy the equations and the boundary conditions. 

To develop the transfer matrix for a lumped mass, a slmlllar pro* 
cedure is followed. The free**body diagram for the lumped mass is shown 
in Fig. 15. Note the assumed form of the response is consistent with 
that of the spring. Since this is a point transfer matrix; 
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Applying Newton's second law of motion and the principle of harmonic 
balance yields the other equations 

R „ L 2 

~ “R 

2 L 2 

- 9 raw x^. 

Note the inclusion of the external excitation !•’ sin((jut) . These can 
now be written in matrix form as 




1 

0 

0 

0 

0 


1 

^1 


2 

— nki) 
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0 

0 

-F 


"1 

X3 

^ - 

0 

0 

1 

0 

0 



A3J 
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0 

-9mw^ 
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.!3 
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0 
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0 
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V. y 


Hence, systems with these types of elements can be modeled by 
appropriate multiplication of transfer matrices and application of 
boundary conditions. However, one problem is quite noticeable. That 
is, the nonlinear terms in the extension column of the spring transfer 
matrix preclude an explicit solution. This problem is overcome in the 
same fashion as finding the roots of the nonlinear polynomials pre- 
viously discussed. 

The solution is initially formulated by ignoring the nonlinear 
terms in the extension column. Once an Initial solution is found, 
this solution is used to predict the magnitude of the nonlinear terms 
in the matrix extension column. This prediction doesn't use the linear 



solution directly, but predicts different values by using a truncated 
Taylor Series. The solution is then reformulated, including the non- 
linear correction terms, and a new solution is found. This process 
is repeated until a desired accuracy is achieved. Again, we will use 
an example to help illustrate the method. 

For the system shown in Fig. 1, the boundary conditions are 

* 0 ; » 0 . 

Hence, the initial solution is found by ignoring the nonlinear 
correction terms in the spring transfer matrix, which results in the 
linear solution 

x^ - F/k(l/(l-Y^) 

2 

2 u) m 

' “ K 

and x^ * 0. Recall that the previous equations developed for the spring 
are 

Xj^+e(3/4 x^^-3/4 x^x^+3/2 x^x^)-^^/k * 0 

Xj +e(-l/4 Xj^^+3/2 Xj^x^+2/k ■ 0 

These equations can be regarded as two simultaneous equations 
for which a solution Is desired. Using the above linear solution as 
an approximate solution to these equations, one can expand these equa- 
tions In a truncated Taylor series yielding two simultaneous linear 
algebraic equations whl<.h can be solved for a more precise solution. 

The resulting solution to these equations can then be used to more 
accurately predict the nonlinear correction terms In the spring transfer 
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matrix. This provides for better convergence of the overall solution. 
The overall solution is then reformulated, using the nonlinear correc- 
tion terms in the spring transfer matrix, and a new solution found. 

This process is repeated until a desired accuracy is achieved. 

In regions where multiple solutions are possible, i.e., as with 
the jump phenomena, one can find these multiple solutions. This is 
accomplished by resolving the problem using a different initial solu- 
tion; i.e., one that differs from the linear solution. This is done 
by modifying the initial state vector in the Initial solution. This 
change is somewhat arbitrary in that it only provides a starting point 
from which to begin the solution process. 

The results obtained (for the system of Fig. 1) with this procedure 
were identical to those obtained using the harmonic-balance principle 
on the differential equation. In obtaining the numerical value of 
the solution, this method required twice the number of iterations as 
the harmonic-balance method. This is because the transfer matrix 
method calculates displacement and internal force Instead of just 
displacement, as in the harmonic-balance method. 

SUMMARY 


This report has outlined what is believed to be very significant 
work in the area of the advancement of transfer matrix methodology. 
Problems with large and small nonllnearltles have been solved with the 
proposed harmonic-balance technique. This technique was subsequently 
Implemented in a transfer matrix solution scheme trhlch allows the solu- 
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tion of nonlinear structural mechanics problems via a transfer matrix 
scheme . 

This work has put the project in the position that provides confi- 
dence that its implementation on a geared rotor analysis can proceed in 
an orderly and successful manner. 


